function [f, F, bigK, W,  bigPi, pi, expsum, sigma_noise, dim, Hyps] = ...
    get_approx_laplace_old(hyps, func, n, n_class, X, y, approxF)



sigma_noise = 1e-7;
dim = length(hyps) / n_class;
if n_class == 0
  disp('n_class == 0');
end
Hyps = reshape(hyps, dim, n_class);
bigK = zeros(n*n_class); K = zeros(n,n,n_class); 
for c = 1:n_class
  K(:,:,c) = feval(func,Hyps(:,c), X, X) + sigma_noise*eye(n);
  bigK(1+(c-1)*n:c*n,1+(c-1)*n:c*n) = K(:,:,c);
end 

f = approxF;
% if nargin <= 6
%   f = alg_3_3(n, n_class, K, y);
% else
%   f = alg_3_3(n, n_class, K, y, approxF);
% end
F = reshape(f, n, n_class);
expsum = repmat(sum(exp(F),2),n_class,1);
%
pi = exp(f)./expsum;
bigPi = zeros(n_class*n,n);
for i = 1:n_class
  bigPi((i-1)*n+1:i*n,:) = diag(pi((i-1)*n+1:i*n));
end
W = diag(pi) - bigPi * bigPi';

return;

